home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 101-125 / disk_121 / dataplot / data drawer / lsq < prev    next >
Text File  |  1992-05-06  |  5KB  |  178 lines
  1.  
  2.  
  3. ' "POLYNOMIAL CURVE FITTER"
  4. ' "FROM WILLIAM G. HOOD, BYTE, JUNE 1987"
  5. ' PRINT "MODIFIED BY DALE HOLT"
  6.  
  7.  
  8.  LN=1000 : LD=11 ' LN=MAX DATA POINTS; LD=HIGHEST DEGREE + 1
  9.  
  10.  N=0: M=0: R2=0: MF=0
  11.  S1=0: S2=0: S3=0: S4=0: P1=0: P2=0
  12.  P3=0: I=0: J=0: J1=0: K=0
  13.  VR=0: MM=0: WT=0: P=0
  14.  DIM X(LN), Y(LN), W(LN), C(LD)
  15.  DIM D1(LD), D2(LD), D3(LD), D4(LD), D5(LD), D6(LD)
  16.  
  17. ' GOSUB MATH 'CALLS THE SUBROUTINE
  18. ' INPUT:
  19. ' N=# OF DATA POINTS
  20. ' X( )=X-COORDINATES OF THE DATA POINTS
  21. ' Y( )=Y-COORDINATES OF THE DATA POINTS
  22. ' W( )=WEIGHTING FACTORS OF THE DATA POINTS
  23. ' M=DEGREE OF POLYNOMIAL
  24. ' MF=0 IF NEW DATA, MF=1 IF OLD DATA BUT HIGHER DEGREE
  25.  
  26. ' OUTPUT:
  27. ' C=ARRAY OF M+1 COEFFICIENTS
  28. ' S2=RESIDUAL VARIANCE
  29. ' R2=COEFFICIENT OF DETERMINATION
  30.  
  31. START:
  32.   CLS
  33.   PRINT "POLYFIT. COPYRIGHT (C) 1986 BY WILLIAM G. HOOD"
  34.   PRINT "MODIFIED IN 1987 BY DALE HOLT"
  35.   PRINT
  36.   PRINT "THIS PROGRAM FINDS THE COEFFICIENTS OF THE NTH DEGREE POLY";
  37.   PRINT "NOMIAL"
  38.   PRINT
  39.   PRINT "Y = C(1)*X^N + C(2)*X^(N-1) + ... + C(N)*X + C(N+1)
  40.   PRINT
  41.   PRINT "THAT FITS A SET OF DATA POINTS IN A LEAST-SQUARES SENSE."
  42.   PRINT"EACH DATA POINT MUST CONSIST OF AN X-VALUE, A Y-VALUE, AND AN
  43.   PRINT"OPTIONAL WEIGHT, SEPARATED BY COMMAS AND TERMINATED BY A";
  44.   PRINT " RETURN."
  45.   PRINT "DATA IS READ FROM A DISK FILE UNTIL THE EOF IS REACHED, ";
  46.   PRINT "OR A"
  47.   PRINT "SPECIFIED NUMBER OF DATA POINTS ARE READ IN FROM THE ";
  48.   PRINT "KEYBOARD."
  49.   PRINT: PRINT "MODIFICATIONS BY DALE HOLT PERMIT SAVING A FILE";
  50.   PRINT " OF POINTS"
  51.   PRINT "GENERATED BY THE POLYNOMIAL WHICH HAS BEEN CALCULATED."
  52.   PRINT "THE ORIGINAL DATA POINTS CAN THEN BE PLOTTED ALONG WITH";
  53.   PRINT " THE CURVES"
  54.   PRINT "GENERATED BY THE LEAST SQUARES FIT FOR A VISUAL CHECK OF FIT
  55.   PRINT
  56.   
  57.   INPUT "NAME OF THE DATA FILE (KYBD: = KEYBOARD)? ",FI$
  58.   INPUT "IS THE DATA WEIGHTED (Y/N)? ",W$
  59.   IF W$<>"Y" AND W$<>"y" THEN W$="N"
  60.   IF FI$<>"KYBD:" THEN L1190
  61. L810:
  62.   PRINT "Keyboard data entry"
  63.   PRINT "How many data points (2 -";LN;")";: INPUT N
  64.   IF N=0 THEN N=LN
  65.   IF N<2 THEN PRINT "INVALID! ";: GOTO L810
  66.   FOR K=1 TO N
  67.     IF W$<>"N" THEN INPUT "X,Y,W? ",X(K),Y(K),W(K): W(K)=ABS(W(K))
  68.     IF W$="N" THEN INPUT "X,Y? ",X(K),Y(K): W(K)=1
  69.   NEXT K
  70. L890:
  71.  PM=LD-1: IF N-1<PM THEN PM=N-1
  72. L900:
  73.   PRINT "Degree of polynomial fit (1-";PM;")? ";
  74.   INPUT "",M : M=ABS(M)
  75.   IF M<1 OR M>PM THEN PRINT "INVALID! ";: GOTO L900
  76.   GOSUB MATH
  77.   PRINT: PRINT "Power Series Coefficients:" 
  78.   FOR J=M+1 TO 1 STEP -1
  79.   PRINT "    ";C(J);TAB(22);"X^";M+1-J
  80.   NEXT J
  81.   PRINT: PRINT "Residual variance: ";S2
  82.    
  83.   PRINT: PRINT "Coefficient of Determination:";R2
  84. L1010:
  85.  INPUT "Make a file of calculated points (Y/N)? ",A$
  86.   IF A$="" THEN A$="N"
  87.   A$=UCASE$(A$)
  88.   IF A$="N" THEN L1130 
  89.   IF A$<>"Y" THEN PRINT "INVALID! ";: GOTO L1010
  90. L1040:
  91.  INPUT "File name? ",F$
  92.  IF F$="" THEN PRINT "INVALID! ";: GOTO L1040
  93.  OPEN F$ FOR OUTPUT AS #5
  94.  FOR I=1 TO N: P=C(1)
  95.    FOR K=1 TO M
  96.      P=P*X(I)+C(K+1)
  97.    NEXT K
  98.    ' PRINT X(I);",";P       ' PRINT TO SCREEN ALSO
  99.    PRINT #5, X(I);",";P
  100.  NEXT I
  101.  CLOSE 5
  102. L1130:
  103.   INPUT "Try another degree? ",A$
  104.   IF A$="" THEN A$="N"
  105.   A$=UCASE$(A$)
  106.   IF A$="N" THEN END
  107.   IF A$="Y" THEN GOTO L900
  108.   PRINT "INVALID! ";: GOTO L1130
  109.  
  110. L1190:
  111.   PRINT "Reading from file: ";FI$
  112.   OPEN FI$ FOR INPUT AS #1
  113.   N=0
  114.   PRINT "X";TAB(15);"Y";TAB(28);"W";:PRINT
  115.   WHILE NOT EOF(1) 
  116.     N=N+1
  117.     IF W$<>"N" THEN INPUT #1,X(N),Y(N),W(N): W(N)=ABS(W(N))
  118.     IF W$="N" THEN INPUT #1,X(N),Y(N): W(N)=1
  119.     PRINT X(N);TAB(14);Y(N);TAB(28);W(N)
  120.   WEND
  121.   CLOSE 1
  122.   PRINT: PRINT N;" points were read from the file."
  123.   GOTO L890
  124.   
  125.   GOTO START
  126.  
  127.  
  128. ' CALCULATION SUBROUTINE
  129. MATH:
  130.  IF MF>0 AND M>MM THEN J1=MM+1: MM=M: GOTO L200
  131.  J1=1: MM=M: S1=0: S2=0: S3=0: S4=0
  132.  FOR I=1 TO N
  133.    WT=W(I)
  134.    S1=S1+WT*X(I): S2=S2+WT: S3=S3+WT*Y(I): S4=S4+WT*Y(I)*Y(I)
  135.  NEXT I
  136.  D4(1)=S1/S2: D5(1)=0: D6(1)=S3/S2: D1(1)=0: D2(1)=1
  137.  VR=S4-S3*D6(1)
  138. L200:
  139.  FOR J=J1 TO MM
  140.    S1=0: S2=0: S3=0: S4=0
  141.    FOR I=1 TO N
  142.      P1=0: P2=1
  143.      FOR K=1 TO J
  144.        P=P2: P2=(X(I)-D4(K))*P2-D5(K)*P1: P1=P
  145.      NEXT K
  146.      WT=W(I): P=WT*P2*P2
  147.      S1=S1+P*X(I): S2=S2+P: S3=S3+WT*P1*P1: S4=S4+WT*Y(I)*P2
  148.    NEXT I
  149.    D4(J+1)=S1/S2: D5(J+1)=S2/S3: D6(J+1)=S4/S2
  150.    D3(1)=-D4(J)*D2(1)-D5(J)*D1(1)
  151.    IF J=>4 THEN 
  152.      FOR K=2 TO J-2
  153.        D3(K)=D2(K-1)-D4(J)*D2(K)-D5(J)*D1(K)
  154.      NEXT K
  155.    END IF
  156.    IF J>2 THEN D3(J-1)=D2(J-2)-D4(J)*D2(J-1)-D5(J)
  157.    IF J>1 THEN D3(J)=D2(J-1)-D4(J)
  158.    FOR K=1 TO J
  159.      D1(K)=D2(K): D2(K)=D3(K): D6(K)=D6(K)+D3(K)*D6(J+1)
  160.    NEXT K
  161.  NEXT J
  162.  FOR J=1 TO M+1
  163.    C(J)=D6(M+2-J)
  164.  NEXT J
  165.  P2=0
  166.  FOR I=1 TO N
  167.    P=C(1)
  168.    FOR J=1 TO M
  169.      P=P*X(I)+C(J+1)
  170.    NEXT J
  171.    P=P-Y(I): P(2)=P2+W(I)*P*P
  172.  NEXT I
  173.  S2=0: IF N>M+1 THEN S2=P2/(N-M-1)
  174.  R2=1: IF VR<>0 THEN R2=1-P2/VR: IF R2<0 THEN R2=0
  175.  RETURN
  176.  
  177.  END
  178.